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Executive  Summary  —  Bubbly  media  play  a  significant  role  in  underwater  acoustics,  medical 
ultrasound  and  in  industrial  systems  where  gas-liquid  flows  are  present.  The  focus  of  our  research 
has  been  to  develop  a  continuum  model  for  bubbly  mixtures  that  can  be  used  to  model  physical 
phenomena  in  these  areas.  The  key  to  the  continuum  model  is  a  nonlinear,  non-equilibrium  equation 
of  state  (EOS)  that  relates  pressure  to  the  mixture  density  and  the  number  density  (number  of 
bubbles  per  unit  volume)  and  their  first  two  material  time  derivatives.  The  derivation  of  the  EOS 
is  presented  here  and  a  number  of  traveling  wave  solutions  obtained  using  this  nonlinear  EOS  are 
discussed.  To  develop  an  accurate  model,  two  important  damping  mechanisms  for  the  medium  had  to 
be  incorporated:  heat  transfer  and  relative  motion  between  the  gas  and  liquid  phases.  To  quantify 
the  importance  of  heat  transfer,  an  analysis  of  single-bubble  radial  oscillations  was  completed  in 
this  work,  and  a  Fade  approximation  for  the  thermal  damping  was  derived  from  the  linearized  gas 
dynamics  equations.  A  second  important  damping  mechanism  arises  from  relative  motion  between 
the  gas  bubbles  and  the  liquid.  The  quantitative  effects  of  relative  motion  on  the  damping  of  waves 
in  bubbly  liquids  has  also  been  examined  and  is  described  here. 
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1  Introduction 


The  study  of  wave  propagation  in  liquids  containing  gas  bubbles  is  an  interesting  multiphase  flow 
problem  because  the  bubbles  themselves  can  behave  in  a  highly  nonlinear  fashion  due  to  their  large 
compressibility.  Bubbly  liquids  occur  in  a  number  of  environmental  and  industrial  settings.  Bubbles 
are  present  near  the  surface  of  the  ocean  in  the  form  of  bubble  clouds  and  plumes.  Air  becomes 
entrained  in  the  ocean  by  breaking  waves  or  the  passing  of  surface  ships.  Bubbles  oscillations  give  rise 
to  underwater  sound  or,  in  the  case  of  surface  ships,  collapse  of  cavitation  bubbles  may  cause  damage 
to  propeller  blades.  In  industry,  the  presence  of  bubbly  mixtures  can  lead  to  inefficient  performance 
in  chemical  reactors  and  boilers,  or  to  safety  problems  for  nuclear  power  plants.  Bubbles  also  play  a 
role  in  applications  of  medical  ultrasound,  such  as  lithotrispy  (breaking  of  kidney  stones  by  sound) 
and  through  their  use  as  contrast  agents  during  ultrasonic  imaging.  The  objectives  of  this  work  were 
to  devise  a  continuum-level  description  of  wave  propagation  and  fluid  flow  in  bubbly  liquids  and  to 
study  nonlinear  waves  (e.g.  shocks)  in  such  media. 

Campbell  and  Pitcher  [1]  carried  out  the  first  qualitative  experiments  on  shock  waves  in  bubbly 
liquids.  They  established  that  the  Mach  number  of  the  wave  must  be  greater  than  unity  for  shocks  to 
exist.  The  Mach  number  was  defined  as  the  speed  at  which  the  shock  propagates  divided  by  the  sound 
speed  of  the  medium  ahead  of  the  shock.  Experiments  performed  by  Noordzij  and  van  Wijngaarden 
[2]  depicted  three  qualitatively  distinct  shock  waveforms  by  stages.  The  first  stage  is  considered  the 
“classical”  shock  profile.  This  profile  can  be  characterized  by  a  steep  rise  in  pressure  from  a  low 
equilibrium  value  followed  by  an  oscillatory  relaxation  region  where  the  waveform  oscillates  about 
a  high  equilibrium  value.  The  oscillations  that  follow  the  initial  discontinuity  are  a  direct  result  of 
the  presence  of  gas  bubbles  in  the  liquid. 

Van  Wijngaarden  [3]  was  the  first  investigator  to  present  a  theoretical  model  describing  bubbly 
liquids  as  a  continuum.  The  derivation  of  the  model  was  primarily  based  on  physical  arguments  and 
so  the  extent  to  which  the  equations  were  valid  was  not  rigorously  defined.  Later,  Caflisch  et  al.  [4] 
started  with  the  equations  describing  the  medium  on  a  microscopic  scale,  and  through  the  use  of 
asymptotic  homogenization,  derived  equations  valid  on  the  macroscopic  scale.  Further  theoretical 
work  wais  done  to  incorporate  additional  physical  effects  into  the  model.  Watanabe  and  Prosperetti 
[5]  included  thermal  conduction,  which  is  considered  the  most  dominant  damping  mechanism  during 
bubble  oscillations.  Most  often  the  bubbles  are  assumed  to  expand  and  contract  adiabatically,  but 
this  is  known  to  be  inaccurate  in  practice.  Sangani  [6]  extended  the  model  to  include  bubble-bubble 
interactions.  Tan  and  Bankoff  [7]  considered  theoretically  the  effects  of  relative  motion  between  the 
two  phases.  Their  work  was  followed  up  a  decade  later  by  Kameda  and  Matsumoto  [8]  and  by  Ishii 
et  al.  [9]  In  these  papers,  the  authors  were  able  to  incorporate  the  effects  of  relative  motion  as  well 
as  some  thermal  effects  into  their  numerical  simulations.  However,  their  results  were  obtained  by 
performing  calculations  on  each  bubble  in  the  flow,  negating  the  advantage  of  considering  a  bubbly 
medium  as  a  continuum.  More  recently,  Kameda  and  Matsumoto  [10]  have  shown  that  there  can  be 
good  agreement  between  theory  and  experiments,  if  the  experiments  are  performed  in  bubbly  liquids 
in  which  the  bubble  distribution  is  spatially  uniform.  There  is  an  extensive  list  of  work  focussing  on 
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shock  propagation  in  bubbly  liquids  which  includes  [ll]-[27]. 

In  order  for  a  continuum  description  of  bubbly  liquids  to  be  valid,  there  needs  to  exist  a  separation 
of  length  scales.  Typically  one  assumes  that  the  distance  between  bubbles  is  much  larger  than  the 
radii  of  the  bubbles,  so  bubble-bubble  interactions  can  be  neglected.  The  size  of  an  averaging 
volume  should  be  large  so  that  it  contains  many  bubbles,  but  small  when  compared  to  the  acoustic 
wavelengths  of  interest.  With  the  latter  assumption,  average  field  quantities  (e.g.  pressure)  can  be 
defined  in  the  averaging  volume,  which  are  spatially  uniform  over  the  volume.  As  a  consequence 
of  this  separation  of  length  scales,  a  non-equilibrium  equation  of  state  (EOS)  can  be  derived.  For 
fiows  where  the  relative  motion  between  the  phases  is  negligible,  this  EOS  is  a  nonlinear  relation 
between  the  pressure  and  the  density  and  its  first  two  material  time  derivatives.  Nigmatulin  [28] 
was  the  first  to  point  out  that  an  EOS  of  this  form  can  be  obtained;  wave  motion  in  such  media  was 
further  investigated  by  Gavrilyuk  [29].  Other  assumptions  for  the  derivation  of  the  EOS  are  liquid 
incompressibility,  monodispersity,  and  absence  of  bubble  breakup  or  coalescence. 

With  the  aim  of  describing  shock  propagation  in  bubbly  liquids  quantitatively,  we  present  several 
ideas  for  theoretical  modeling  and  discuss  preliminary  results.  We  first  consider  the  issue  of  thermal 
damping  of  a  single  bubble  in  an  infinite  liquid.  Insightful  work  investigating  thermal  damping 
in  single  bubble  oscillations  can  be  found  in  papers  by  Devin  [30]  and  Prosperetti  [31].  Solutions 
for  the  equation  of  motion  for  an  isothermal  or  an  adiabatic  bubble  are  straightforward  but  not 
very  accurate  in  modeling  real  bubbles.  By  accurately  describing  heat  conduction  for  a  linearly 
oscillating  bubble,  we  develop  an  effective  thermal  viscosity  which  will  contribute  to  the  damping 
of  sound  waves  in  bubbly  liquids.  Values  of  this  viscosity  are  found  to  be  only  a  function  of  the 
equilibrium  radius  of  the  bubble  and  can  be  written  in  the  form  of  a  Fade  approximate. 

In  the  next  phase  of  our  analysis  we  derive  the  non-equilibrium  EOS  for  a  monodisperse  bubbly 
liquid  incorporating  the  relative  motion  between  the  gas  phase  and  the  liquid  phase.  The  importance 
of  the  EOS  is  that  it  contains  all  the  physics  needed  to  characterize  the  nonlinearity  of  the  medium 
in  a  continuum  model.  This  approach  is  meritorious  both  theoretically  and  computationally  for  its 
simplicity.  We  will  not  have  to  keep  track  of  each  individual  bubble  in  order  to  calculate  nonlinear 
wave  propagation  through  the  bubbly  medium.  Including  relative  motion  into  the  model  leads  to  a 
greater  understanding  of  the  damping  mechanisms,  and  how  they  interact  with  one  another,  within 
bubbly  liquids.  Also  presented  herein  are  some  results  which  are  obtained  through  this  approach.  By 
seeking  a  traveling  wave  solution  using  the  EOS  in  conjunction  with  the  fully  nonlinear  continuity 
and  Euler  equations,  a  phase-plane  description  of  the  system  can  be  obtained.  The  phase  portrait 
provides  insights  into  the  qualitative  behavior  of  shocks  in  bubbly  liquids. 

2  Thermal  Damping  of  a  Single  Bubble:  Fade  Approxima¬ 
tion 

The  derivation  for  the  thermal  damping  of  a  single  oscillating  bubble  begins  with  the  Rayleigh-Plesset 
equation.  The  Rayleigh-Plesset  equation  describes  the  nonlinear  motion  of  a  radially  oscillating 
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bubble  in  an  infinite  liquid.  It  is  a  second  order  ordinary  differential  equation  (ODE)  for  the  radius 
R{t)  of  the  bubble  as  a  function  of  time  t: 

pt[RR+lR^]  +  4,ie^  =  Pgw-P°°(t).  (1) 

Here  pt  and  pt  are  the  density  and  viscosity  of  the  liquid,  is  the  forcing  pressure  far  from  the 
bubble  center  and,  for  this  derivation,  it  is  considered  to  be  known.  Pgw  is  the  pressure  inside  the 
bubble  evaluated  at  the  bubble  wall.  In  order  to  solve  Eqn.  (1),  the  pressure  inside  the  bubble  must 
be  determined.  For  this  purpose  the  equations  of  gas  dynamics  within  the  bubble  need  to  be  solved. 
They  are: 


=  0 

(2) 

=  0 

(3) 

DT  DP 

Dt  Dt 

=  V-[«:VT] 

(4) 

P 

=  pUT 

(5) 

representing  the  conservations  of  mass,  momentum,  and  energy  and  the  equation  of  state  for  an 
ideal  gas,  respectively.  In  the  above  equations  p,  P,  v  and  T  are  the  density,  pressure,  velocity  and 
temperature  of  the  gas.  Whereas  «  is  the  thermal  conductivity,  Cp  is  the  specific  heat  capacity  of 
the  gas  at  constant  pressure,  and  Tl  is  the  universal  gas  constant.  Note  that  the  viscosity  of  the  gas 
is  neglected  in  this  analysis  whereas  thermal  conduction  is  included.  We  now  assume  a  perturbation 
expansion  of  Eqns.  (l)-(5)  in  the  following  form 

p  =  po  -h  epi  +  . . . 

P  =  Pq  cPi  “h  . . . 

V  =  0  +  evi  +  . . . 

T  =  ro-heTi  +  ... 

R  =  Rq  +  cRi  +  . . . 

The  0{e)  equations  in  the  radial  direction  are 


dt  dr  ^ 

’^Ul) 

=  0 

(6) 

dvi 

dPi 

=  0 

(7) 

dr 

dTi 

dPi 

K  d  j 

( 2dTi\ 

(8) 

dt 

dt 

J.2  Qj.  1 

I  dr) 

Pi 

=  P'iToPi  +  PoPi) 

(9) 
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with  the  Rayleigh-Plesset  written  as 


piR^Ri  +  ”  -Pilr=i  ~  —P'^{t) 


We  scale  the  variables  as  follows 


t*  =  U}„t 


"  To  ‘  RoiOc 

The  time  scale  is  based  on  the  period  of  oscillation  for  a  single  adiabatic  bubble  whose  natural 
frequency  can  be  written  as 

°  ^  PtRl' 

Rewriting  the  0(e)  equations  in  dimensionless  form  (dropping  the  *’s)  yields 


^  1  ^  _  Q 

dt  Ml  dr 


dt  ^  j  dt  Pe  dr  \  dr 


Pi  =  (Pi+Ti)  (14) 

In  Eqns.  (11)“(14)  we  have  introduced  two  key  nondimensional  parameters,  the  Mach  number  and 
the  Peclet  number: 

Mb  =  u;oR<ii[^ 

V  ^0  f^/pcp 

The  Peclet  number  Pe  is  a  measure  of  convective  to  diffusive  effects.  If  the  Peclet  number  is  small, 
there  is  a  thick  thermal  boundary  layer  in  the  bubble  and  the  bubble  can  be  considered  isothermal. 
For  large  Peclet  numbers,  on  the  other  hand,  the  bubble  has  a  thin  thermal  boundary  layer  and  can 
be  considered  adiabatic.  The  second  dimensionless  parameter  is  the  acoustic  Mach  number,  M^, 
inside  the  bubble.  For  an  ideal  gas  inside  an  oscillating  bubble,  the  acoustic  Mach  number  is  small. 
Therefore,  the  r-momentum  equation  reduces  to 


Pi  =  Pi{t) 
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We  add  the  continuity  and  energy  equations  and  use  the  ideal  gas  equation  of  state  Eqn.  (14)  to 


=  (15) 

T  dt  ■'  Per’ Sr  V 

Multiply  Eqn.  (15)  by  and  integrating  from  0  to  1  yields 

r^dPi\  2  |i  1  2dTi^  (.f.s 

Eqn.  (16)  gives  an  expression  for  time  derivative  of  the  pressure  in  terms  of  the  heat  flux  at  the 


surface: 


dt  ^  ^  Pe  dr 


Equation  (17)  can  now  be  used  in  the  energy  equation  (13),  resulting  in  an  equation  for  the  tern- 
perature  field  only: 

at  Pe  dr  J  Pe  r"^  dr  \  dr  J  ’ 

Since  the  gas  is  in  thermal  equilibrium  with  the  liquid  and  the  temperature  of  the  liquid  is  a  constant, 
To,  we  have  the  following  boundary  for  Ti  for  the  equation  above 

Ti  =  0  at  r  =  1 .  (19) 

We  now  seek  solutions  to  (18)  of  constant  frequency 

Pi  =  3i[Pe’”‘] 


Ti  =  5R[T(r)e‘' 


Pi  =  9?[Pe’”‘] 

Eqn.  (18)  becomes 

inf  +  3(-,  -  1)  [iilfl  -  if  (1)]  =  ')  (20) 

Since  the  second  term  in  Eqn.  (20)  is  a  constant,  we  define  a  new  variable,  0,  which  takes  this  into 
account 

0  =  T  +  3(7-1)  •  (21) 

Now  we  have  a  much  simpler  ODE  to  solve,  namely 

‘ne(’')  =  ^^|(--“0').  (22) 
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with  the  boundary  conditions 


0(1)  =  3(7-1) 


0(0)  is  finite 

For  differential  equations  of  this  form,  rewriting  the  dependent  variable,  0,  as  the  ratio  can 
lead  to  the  simplified  ODE 


2i\Q{r)=Q”{r) 


in  which 


with  new  boundary  conditions 


Q(l)  =  3(7-1)  R+^T'{1) 


Q(0)  =  0 


Eqn,  (23)  has  the  general  solution 


Q{r)  =  A  sinh[(l  +  i)Xr]  +  B  cosh[(l  4-  i)Xr] . 


Applying  the  boundary  conditions  gives 


Qir)  -  3(7  1)  ^  (1)J  sinh[(l  +  i)X]  ' 


We  now  use  Q  =  Q/r  and  Eqn.  (21)  to  write 


iJr+  *  /'sinh[(l  +  i)Ar] 

^  ^^JVsinh[(l  +  i)A]  ^ 


f  =  3(7- 


An  expression  for  T''(l)  is  found  by  taking  the  derivative  of  Eqn.  (25)  at  r=l. 

3(7  -  1)  [(1  +  i)Acoth[(l  +  i)A]  -  1]  B 
T  (1)  =  : - 3(7-ib  r..  7" — ..r.,  .  .s,, 


1  -  [(1  +  *)A  coth[(l  +  f)A]  -  1]  . 

Rrom  Eqns.  (18)  and  (26)  we  have  an  expression  for  the  pressure  inside  the  gas  bubble  at  the  wall. 

iUP  =  +  I^T'(l) 


P-  37i?  *2A2  [l-i5^$(A)J 


$(A)  =  (1  +  i)Acoth[(l  +  i)A]  -  1 . 


In  the  limit  of  large  Pe,  A  oo  and  $  simplifies  to 

$(A)  =  (l  +  i)A-l. 

Given  the  pressure  inside  the  bubble  as  a  function  of  frequency,  it  is  possible  to  write  down  the 
response  of  an  oscillating  bubble  from  the  Rayleigh-Plesset  equation  (1): 

+  i2CflR  +  3(^  _  +  jHPe^  ^  ^ 

This  is  the  exact  linear  solution  for  an  oscillating  bubble.  For  the  unforced  case,  the  free  response 
for  a  bubble  has  a  frequency  which  is  a  root  of: 

+  ^2(^n  +  ^  =  0  (29) 

The  thermal  effects  for  an  oscillating  bubble  are  contained  in  the  last  term  in  Eqn.  (29),  which  comes 
from  the  pressure  inside  the  bubble  and  can  be  rewritten  more  clearly  as 

The  real  part  of  the  pressure  term  provides  the  natural  frequency  of  the  bubble.  This  natural 
frequency  includes  a  shift  that  takes  into  account  thermal  effects.  One  half  the  imaginary  part  of 
the  pressure  term  is  the  contribution  to  damping  by  heat  conduction.  For  comparison,  the  damping 
coefficient  due  to  liquid  viscosity,  near  resonance  is  1.27  x  10“®  for  a  100  micron  bubble  in  water 
(10°G,  1  atm).  The  damping  coefficient  due  to  thermal  effects  is  42.7  x  10“^. 

The  solution  for  fi  from  Eqn.  (29)  can  be  expanded  for  both  large  and  small  Pe.  For  large  Pe 
define  a  small  parameter,  e,  as 


and  take  fl  to  be 


n(Pe  -}•  oo)  =  =  Go  +  efii  +  + 


After  taking  the  expansion  described  in  Eqn.  (30)  and  substituting  it  into  Eqn.  (29),  we  obtain  the 
following  set  of  algebraic  equations: 


0(1): 


1  -  fig  =  0 


0(£): 


0(e2):  T  (2fio)3/2 

Note  that  we  scale  C  with  the  perturbation  parameter,  e,  so  that  ^  =  ^e.  Taking  the  leading  order 
solution  Q.Q  =  lj  the  first  two  corrections  of  are 

«■  ■ 

fi2  =  -g(i  +  l)  (3(l  +  ^)(37-l)(T-l)-3^/2(7-l)C  +  2(l  +  ^)^) 


-2fio(fii-tC)-^^'^  -'^=0 


-fij  —  2  J7o  H"  ^  2  0,1  C  ~  ^ 


.3(37-2)(7-1)  ,  3(7-1)(1-0»i 
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For  small  Pe  the  complete  expansion  of  fl  is  thus: 


n(Pe  —^0)  =  fl®  =  Wo  +  e  W2  +  •  •  • 


(31) 


in  which  the  small  parameter,  e  is  now 


e  =  Pe 


Eqn.  (29)  with  Eqn.  (31)  yields  the  algebraic  equations 

0(1): 


i-w2  =  0 

7 


0(e): 


0(6^) 


(37  +  7)(7-l)cug 


2  u)\  —  0 


.  Wi  (7  -  1) 


—  2  Wo  W2  +  i  — — ~  +  *  2  ^  Wl  —  Wj 


15757®  — "  15^2 

Taking  the  leading  order  solution  wo  =  ^,  the  first  two  corrections  to  are 


Wi 


W2 


_  (7-1)  (7 +  7)  C(7-l) 

25207^/2  30  73/2  2 


Since  we  are  interested  in  the  damping  due  to  thermal  effects  and  not  viscosity,  we  will  now  set 
^  =  0.  The  first  two  corrections  for  Eqn.  (29)  are  known  in  terms  of  Pe  for  both  large  and  small  Pe 
numbers. 


1  1 
SS  flo  +  fil  /=—  +  ^ 

VPe 

«  Wo  +  wi  Pe  +  W2  Pe® 

These  expansions  can  be  combined  in  a  two-point  Pade  approximation  of  the  form 

<-iyn  ^  \  oo  +  oi  Pe^^^  +  02  Pe 

n(Pe)  «  n(Pe)  = - TT? - 

^  ^  ^  ^  l-l-6iPe^/®-t-62Pe 

With  the  Padd  expression  above,  only  the  first  two  terms  from  the  expansion  in  Eqns.  (32)-(33)  can 
be  matched.  Expanding  Eqn.  (34)  for  both  large  and  small  Peclet  numbers  gives 


(32) 


(33) 


(34) 


For  Pe  0 
n(Pe) 
For  Pe  00 
n(Pe) 


oo  +  (oi  -  oo  61)  -I-  (02  -  oi  61  -f  Oo  {h\  -  62))  Pe 

1 


02  /  02^1 

52  V  51  b^J 
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The  coefficients  from  the  Fade  expression  can  be  found  by  matching  them  to  the  coefficients  from 
the  solution  expansions  (32,33). 


(fl  -  (fi  +  n*)  =  0  (35) 

Here  is  the  complex  conjugate  of  Cl.  We  denote  the  damping  due  to  thermal  effects  by  Qh-  From 
(35)  we  can  write 

Fig.  1  is  a  plot  of  the  magnitude  of  the  response  of  a  bubble  being  forced  periodically.  The  plot  is 
for  a  Pe=  25,  corresponding  roughly  to  a  bubble  with  a  radius  of  25  /xm.  Shown  in  the  plot  are  four 
curves  representing  solutions  obtained  from  the  “exact”  equation  (28),  from  the  Fade  approximation, 
and  from  the  adiabatic  and  the  isothermal  limits  of  Eqn.  (28).  For  a  Pe=  25  the  plot  shows  that 
the  adiabatic  solution  underdamps  the  response  while  the  isothermal  solution  overdamps.  Our  Pade 
approximation  provides  more  accurate  results  for  the  range  of  Pe  numbers  that  we  are  interested  in 
corresponding  to  bubble  radii  between  5-100  ^m. 

3  The  Non-Equilibrium  Equation  of  State  (no  relative  mo¬ 
tion) 

We  now  consider  a  dilute  mixture  of  identical  spherical  gas  bubbles,  having  radii  i?(x,i),  which  are 
suspended  in  a  liquid  medium  of  density  pi  and  viscosity  p.  When  a  separation  of  length  scales 
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15. 


Pe  =  25 


Figure  1: 


holds  such  that  the  bubble  radii  are  much  smaller  than  the  typical  distance  between  bubbles,  which 
is  in  turn  smaller  than  the  size  of  an  averaging  volume  within  which  a  large  number  of  bubbles 
can  be  found,  a  non-equilibrium  equation  of  state  which  relates  pressure  and  density  in  the  mixture 
can  be  obtained.  This  also  requires  the  typical  wavelength  of  sound  waves  in  the  mixture  to  be 
large  compared  to  the  size  of  the  averaging  volume,  so  that,  to  a  good  approximation,  the  acoustic 
pressure  can  be  considered  to  be  spatially  uniform  on  the  scale  of  the  averaging  volume.  Under 
these  approximations,  the  equation  of  motion  for  the  radial  oscillations  of  each  of  the  identical 
noninteracting  bubbles  in  the  averaging  volume  is  simply  the  Rayleigh-Plesset  equation  [32] 

Pe[RR  +  =  -P(x,  t) .  (36) 

This  is  a  nonlinear  second  order  ODE  for  the  bubble  radius  as  a  function  of  time.  Here,  each  overdot 
represents  a  time-derivative,  Po  and  Ro  are  the  equilibrium  values  of  the  pressure  in  the  mixture  (and 
in  the  bubbles)  and  the  bubble  radii,  respectively,  and  P{x,  t)  is  the  mean  pressure  in  the  mixture  at 
the  position  of  the  averaging  volume.  Effects  of  surface  tension  have  been  neglected,  while  viscous 
damping  of  bubble  pulsations  has  been  taken  into  account.  In  (36),  7  denotes  the  poly  tropic  index 
of  the  gas  inside  the  bubble;  its  value  ranges  from  unity,  for  isothermal  bubble  oscillations,  to  the 
ratio  of  constant-pressure  to  constant-volume  heat  capacities  of  the  gas  for  adiabatic  oscillations. 
Generally  in  the  Rayleigh-Plesset  equation,  the  forcing  pressure  on  the  right-hand  side  is  regarded 
as  the  pressure  far  from  the  individual  bubble;  due  to  the  assumed  separation  of  scales,  this  pressure 
is  approximately  the  same  as  the  mean  pressure  in  the  averaging  volume  surrounding  the  bubble. 

The  volume  fraction  <j)  of  bubbles  in  the  mixture  can  be  related  to  their  number  density  n  (number 
of  bubbles  per  unit  volume)  by 

<j)  =  .  (37) 

o 

In  terms  of  <l>  the  density  of  the  bubbly  mixture  is  given  by 

p  =  pt{\  -<!>)+  Pg<l>  »  PtO-  -  <A)  • 


12 


In  the  last  approximation,  we  have  neglected  the  contribution  of  the  gas  phase  (density  pg)  to  the 
mass  density  of  a  dilute  bubbly  liquid.  In  the  absence  of  bubble  breakup  and  coalescence,  the  number 
of  bubbles  per  unit  mass  of  mixture  remains  constant  with  time,  i.e. 

(39) 


n  .  . 

—  =  constant  =  —  . 
P  Po 


The  subscript  ‘o’  refers  to  the  equilibrium  state.  Eqs.  (37)  through  (39)  can  now  be  combined  to 
yield  a  one-to-one  relationship  between  the  mixture  density  p  in  the  averaging  volume  and  the  radii 
of  bubbles  in  that  volume: 

f  =  l-<Ao(|-)^f  •  (40) 

pi  po 

Upon  solving  Eq.  (40)  for  J?  as  a  function  of  p  and  substituting  the  result  into  the  Rayleigh-Plesset 
equation  (36),  one  obtains  a  relationship  of  the  general  form 


P  =  Po 


_Pt^o_P_V  ^  ^PPt  Dp 

Po{pi-p)\  ip{Pi-p)Dt 


peRliPol4>o)y^ 

3p^{l/p-l/pty/^ 


DP 


Pi 


mi 


Eq.  (41)  relates  the  mixture  pressure  P  to  the  density  p  and  its  first  two  time-derivatives.  Pro¬ 
vided  that  translational  motion  of  bubbles  relative  to  the  surrounding  liquid  is  negligible,  the  time- 
derivatives  in  (41)  are  interpreted  as  material  derivatives  when  the  spatial  variation  of  P  and  p  on  the 
macroscale  is  taken  into  account.  A  more  in  depth  derivation  of  this  fully  nonlinear,  non-equilibrium 
equation  of  state  can  be  found  in  the  paper  by  Nadim,  Goldman  and  Barbone  [33]. 


3.1  Traveling  Wave  Solution 

In  a  continuum  description  of  any  fluid  mixture  in  which  direct  viscous  dissipation  effects  are  neg¬ 
ligible,  the  laws  of  conservation  of  mass  and  linear  momentum  in  the  continuum  take  the  standard 
forms: 


dp  dp  du 
dt  ^  ^  dx  ^  dx 


,du  du.  dp 


=  0, 
0. 


(42) 

(43) 


Here,  p  and  p  refer  to  the  mean  density  and  pressure  in  the  mixture  and  u  denotes  the  x-component 
of  the  translational  velocity.  Here,  only  one-dimensional  motions  along  spatial  coordinate  x  are 
considered  and  all  three  field  variables  are  dependent  upon  {x,t)  where  t  represents  time. 

One  can  seek  a  traveling  wave  solution  of  the  above  equations  by  defining  the  traveling-wave 
coordinate  p  as 

7]  =  X  —  Ut,  (44) 

U  is  the  propagation  speed  of  the  waveform.  A  coordinate  transformation  from  (x^t)  to  r}  reduces 
(42)  and  (43)  to  the  coupled  system  of  ODEs 

{u-U)p^  +  pu^  =  0,  (45) 

p(u-I/)u'-hF'  =  0.  (46) 
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A  prime  represents  a  total  derivative  with  respect  to  Tj.  Eqs.  (45)  and  (46)  can  each  be  integrated 
once  to  yield 


piu-U)  =  Cl,  (47) 

Ciu  +  P  =  C2.  (48) 


Cl  and  C2  are  constants  of  integration. 

We  suppose  that  at  large  positive  x,  i.e.,  as  rj  tends  to  infinity,  the  velocity  u  in  the  medium  is 
zero  and  the  pressure  and  density  attain  constant  equilibrium  values  po  and  po,  respectively.  The 
integration  constants  Ci  and  C2  can  be  thus  evaluated  based  upon  conditions  at  infinity  to  be 
Cl  =  —poU  and  C2  =  Po-  When  these  values  are  substituted  into  (47)  and  (48)  and  the  variable  w  is 
eliminated  from  the  two  equations,  the  result  is  a  single  nonlinear  algebraic  equation  which  relates 
the  pressure  and  density  profiles,  P(»?)  and  p{t))‘. 

P  =  Po  +  PoU^i^^).  (49) 

This  relationship  is  exact  for  one-dimensional  traveling  waves  in  any  continuum  which  is  described 
by  the  standard  continuity  and  Euler  equations.  For  the  bubbly  liquids  which  are  the  subject  of 
this  study,  one  can  find  a  further  relationship,  in  the  form  of  a  non-equilibrium  equation  of  state, 
between  pressure  and  density  (and  its  first  two  derivatives)  in  the  mixture.  The  latter  can  be 
combined  with  the  nonlinear  algebraic  equation  (49)  to  reduce  the  problem  to  a  single  nonlinear 
ODE  that  is  amenable  to  phase-plane  analysis. 


In  the  traveling-wave  coordinate  7/,  the  equation  of  state  (41)  takes  the  form 

Pt4>oP  V 


P  =  Po 


AppipoU 


{Po{pt-p)\  ip‘^{pi-p) 


,  PtRl{Pol<t>o?f^ 


( 


9l 

6p(p^-p) 


(50) 


When  this  equation  is  combined  with  the  nonlinear  algebraic  relation  (49)  resulting  from  the  inte¬ 
gration  of  continuity  and  Euler  equations,  a  second-order  nonlinear  ODE  is  obtained  for  the  density 
profile  p(7j).  The  latter  ODE  is  most  simply  given  in  its  dimensionless  form.  For  this  purpose,  we 
define  dimensionless  variables 


P*  =  P/Po 
P*  =  P/Po 

77*  = 

p\  =  P^/Po  =  1/(1  -  0o) 


where  quantities  with  superscript  are  dimensionless.  If  this  superscript  is  now  dropped  for  nota- 
tional  simplicity,  the  resulting  dimensionless  nonlinear  ODE  for  the  density  profile  of  the  traveling 


wave  is 


1+^M2 

00 


Pi  <l>o  p  _  yCpt 

.pi-p\  p^{pt-p)^ 
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(51) 


+ 


1/3 

IPl 


<^T  p”/®  {.pt  -  pyi^  i  \^p'^p<-  -  p) 


p"+ 


(- 


pt 


The  terms  that  appear  on  the  left-hand  side  of  Eqn.  (51)  come  from  the  integrated  form  of  the  Euler 
equations  and  the  terms  on  the  right-hand  side  are  from  the  EOS.  Two  dimensionless  parameters, 
^  and  M,  fully  characterize  this  system.  These  are  defined  by 

4/Li 

^  ~  PtRl^o' 


M  =  U/co. 

The  parameter  C  represents  a  dimensionless  damping  constant  while  the  Mach  number  M  is  the 
ratio  of  the  speed  of  the  traveling  wave  to  the  low-frequency  sound  speed  in  the  bubbly  liquid.  In 
Eq.  (51),  a  prime  denotes  a  derivative  with  respect  to  the  dimensionless  traveling-wave  coordinate 

V- 


3.2  Phase-Plane  Description  of  Traveling  Waves 


In  order  to  explore  the  qualitative  features  of  the  solutions  of  Eq.  (51),  it  is  convenient  to  consider  its 
phase-plane.  The  use  of  a  phase-plane  description  in  qualitatively  describing  traveling  shock  waves 
and  the  nature  of  the  fixed  poiiits  can  also  be  found  in  Tan  and  Bankoff  [7^;  a  similar  phase-plane 
analysis  of  traveling  waves  in  an  Euler-Poisson  model  of  a  plasma  is  presented  by  Cordier  et  al.  [34]. 

To  examine  the  phase  plane,  the  second-order  nonlinear  ODE  is  written  as  a  system  of  two 
first-order  equations  by  introducing  the  new  variable  g  to  denote  p'.  This  results  in  the  coupled 
system 


9, 

A{p) 

l  +  1 

00 

fp-i- 

\  (Pl<t>oP'^ 

1  p  . 

)  \Pt-p/ 

'  p^  {pt  -  p) 

(1- 

Pt  \ 

1 

_  P^  \  a- 
V/>  ^p{pi-p))^ 


(52) 


(53) 


Where, 


A{p)  = 


The  phase-plane  displays  all  possible  solutions  of  this  system  of  equations,  with  p  on  the  horizontal 
axis  and  g  (i.e.  p')  on  the  vertical  axis.  To  find  the  fixed  points  of  the  system,  we  set  both  p'  and  g' 
to  zero.  The  following  equation  describes  the  locations  of  the  fixed  points. 


1  +  -^ 
(Po 


/p-1^  _  \pe<l>oPy  _Q 
\  P  )  [pe-p\ 


(54) 
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Since  pi  =  l  satisfies  Eqn.  (54)  it  is  a  fixed  point.  Note  that  pi  -  1/(1  -  ^o)-  The  other  fixed  point 
can  be  found  graphically  by  plotting  Eqn.  (54).  Fig.  (2)  is  a  plot  showing  the  two  cases  possible, 
M  >  1  and  M  <  1.  For  M  >  1  the  solid  curve  crosses  the  p-axis  at  a  value  of  p  higher  than  unity. 
This  is  the  location  of  the  second  fixed  point.  For  M  <  1  there  is  no  second  fixed  point,  indicated 
by  the  fact  that  the  dashed  curve  never  crosses  the  p-axis. 


Figure  2: 


Figure  3: 


Fig.  3  provides  a  typical  phase  portrait  of  the  system  for  the  case  when  Mach  number  exceeds 
unity  —  in  this  case  M  =  1.5.  As  seen  in  this  figure,  the  fixed  point  at  pi  =  1  is  a  saddle  point,  with 
trajectories  in  the  first  and  third  quadrants  moving  away  from  this  point  and  those  in  the  second 
and  fourth  quadrants  coming  towards  it,  for  increasing  77.  The  other  fixed  point  is  a  spiral  node 
(when  C  =  0.1).  The  location  of  the  second  fixed  point,  p2,  is  the  second  root  of  Eqn.  (54).  There 
is  a  single  trajectory  which  connects  the  two  fixed  points,  identified  in  Fig.  3  as  the  shock  solution. 
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This  trajectory  spirals  away  from  the  fixed  point  at  p2j  eventually  tends  toward  the  saddle  point 
at  p  =  1,  for  increasing  rj.  The  solution  profile  p(ry)  corresponding  to  this  trajectory  is  displayed  in 
Fig,  4. 


C  =  o.i  M  =  1.5 


Figure  4: 


This  solution  represents  a  density  disturbance  of  steady  shape  that  propagates  to  the  right  at 
dimensionless  traveling-wave  speed  M  =  J7/co  =  1.5,  The  dimensionless  density  at  large  positive 
7}  approaches  unity  (i.e.,  the  density  approaches  its  equilibrium  value  po  far  ahead  of  the  shock,  as 
required  by  the  original  boundary  conditions  at  large  positive  x).  At  large  negative  rj,  far  behind  the 
shock,  the  density  settles  down  to  the  higher  value  p2  after  a  series  of  oscillations  which  are  associated 
with  the  spiral  trajectory  near  the  node  in  the  phase  plane.  The  fact  that  these  oscillations  exist  at 
all  is  related  to  the  presence  of  the  second  derivative  term  in  ODE  (51),  which  is  in  turn  associated 
with  the  second  derivative  term  in  the  non-equilibrium  equation  of  state,  (41),  and  the  Rayleigh- 
Plesset  equation  (36).  Thus,  the  oscillations  observed  behind  shock  waves  in  bubbly  liquids  are 
simply  a  result  of  the  volume  oscillations  of  gas  bubbles  in  the  medium  when  they  are  subjected  to 
a  somewhat  sudden  increase  in  pressure,  as  the  pressure  wave  goes  by. 

It  should  be  noted  that  the  direction  of  the  arrows  in  Fig.  3  corresponds  to  increasing  rj.  The 
trajectories  in  the  phase-plane  can  also  be  interpreted  as  functions  of  time  but  with  their  directions 
reversed.  For  instance,  as  a  function  of  time,  at  a  fixed  spatial  position  in  space  initially  far  ahead 
of  the  shock,  the  density  would  start  near  its  equilibrium  value  pi  =  l  (the  saddle  fixed  point)  and, 
as  time  progresses,  move  along  the  shock  trajectory  (in  the  sense  opposite  to  those  indicated  by  the 
arrows)  and  end  up  at  the  higher  density  p2  for  large  positive  times. 

4  Equation  of  State  (with  relative  motion) 

Relative  motion  has  been  shown  in  papers  by  Kameda  and  Matsumoto  [10]  and  Ishii  et  al.  [9] 
to  be  important  in  determining  the  exact  waveform  of  shocks  in  bubbly  liquids.  It  is  possible  to 
derive  a  new  EOS  that  allows  for  relative  motion  between  the  gas  phase  (bubbles)  and  the  liquid 
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phase.  When  there  was  no  relative  motion  between  the  phases,  it  was  possible  to  write  the  bubble 
radius,  i?,  as  a  function  only  of  the  mixture  density [33].  If  we  allow  for  relative  motion,  the  exact 
functionality  of  R  becomes  more  complicated.  The  number  of  bubbles  that  are  in  a  given  volume 
must  simultaneously  be  tracked.  This  leads  to  a  new  relation  for  the  bubble  radius. 


R{p,n) 


Pt-P 

^pen' 


(55) 


Here  n  is  called  the  number  density  £md  is  defined  as  the  number  of  bubbles  per  unit  volume.  Both  p 
and  n  are  field  variables  and  depend  on  time.  If  Eqn.  (55)  is  substituted  into  Eqn.  (36),  the  resulting 
expression  for  the  average  pressure  is 


P(x,  t) 


np  +  (pt  -  p)n 
6(^/9f)tnt(p<-p)5 


_ 5np _ 

%^pt)^ni{pi-  p)  3 


ll{pt-  p)^n^ 
18(^p^)  3  nt 


18(fp^)tnf(p^-p)t 


+  4/i 


np  +  (pt  -  p)n 
Zn{pi  -  p) 


\no{pi-p)J 


(56) 


where  the  overdots  represent  the  convective  derivative  with  respect  to  the  gas  phase  velocity.  In 
Eqn.  (56),  the  first  term  in  the  bracket  is  from  inertial  forces,  the  second  term  is  due  to  viscous 
stresses  on  the  bubble  interface  and  the  last  term  is  from  the  ideal  gas  law  for  adiabatic  bubbles. 
We  now  need  equations  describing  how  the  number  density,  n,  and  the  gas  phase  velocity,  v,  change 
with  space  and  time.  For  n  we  can  use  a  continuity  equation  of  number  density  as  long  as  we  assume 
that  bubbles  cannot  be  created  or  destroyed  during  our  time  of  interest.  For  v  we  will  use  the  gas 
phase  momentum  equation  which  reduces  to  a  force  balance  on  a  bubble  if  the  inertia  of  the  gas  is 
ignored.  These  equations  are 

^  +  V  •  (nv)  =  0 .  (57) 


0 


1  DgR 
R  Dt 


(v  -  u)  -  -pi(j) 


(£^ 

V  Dt 


Dt^\ 
Dt  ) 


—  127r/i£i?n(v  -  u)  -f  picj) 


(58) 


The  forces  in  Eqn.  (58)  were  multiplied  by  n,  to  obtain  the  total  force  per  unit  volume.  The  first 
term  in  Eqn.  (58)  is  the  force  resulting  from  changes  in  the  bubble  radius  as  it  translates,  and  the 
last  three  terms  in  Eqn.  (58)  are  due  to  the  added  mass,  the  drag  and  buoyancy,  respectively.  Note 
that  in  the  force  balance  equation  there  exist  two  different  material  time  derivatives,  following  either 
the  bubble  velocity,  or  the  mixture  velocity: 
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Dt 
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If  in  Eqn.  (58)  buoyancy  is  neglected  and  Eqn.  (55)  is  used  to  replace  R,  then  we  will  have  a 
equation  describing  the  velocity  of  the  gas  phase  in  terms  of  the  mixture  density,  mixture  velocity 
and  the  number  density: 


Dt 


=  3^  +  (v-u) 


^  ^gP  ^ 
pe-p  Dt 


1  DgTl 
n  Dt 


-  18i/f 


‘inptn 
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(59) 


4.1  Traveling  Wave  Solution 

In  a  manner  similar  to  the  case  of  no  relative  motion,  we  seek  a  traveling  wave  solution  to  the  system 
of  equations  describing  bubbly  liquids  with  relative  motion.  As  before,  the  governing  equations  are 
recast  in  terms  of  the  traveling  wave  coordinate  77,  where  7]  =  x  —  Ut.  The  conservation  equations 
are: 


{u-U)p'+pu^  =  0, 

(60) 

{v  —  ?7)n'  +  nv*  =  0 , 

(61) 

p(u-[/)u'+P'  =  0, 

(62) 

After  integrating  these  equation  we  can  assume  for  the  boundary  condition  at  77  ^  00  (x  00)  that 
the  velocities  are  zero  and  the  pressure,  mixture  density  and  number  density  go  to  their  equilibrium 
values  of  Po,  Po  and  n^,  respectively.  From  these  results  we  have  equations  relating  the  velocities 
and  pressure  to  the  density  and  number  density: 


u 

p 

(63) 
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II 

0 

(64) 

p 

(65) 

Equations  (63-65)  can  be  combine  with  the  traveling  wave  versions  of  Eqns.  (59)  and  (56)  to  obtain 
two  ODEs  for  p(77)  and  7^(77).  For  brevity  we  only  provide  the  system  of  ODEs  in  matrix  form  as: 
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^  ^  n^p(pt-p) 

9 

.  0 

0  n 

1 

,  f2{p,n,g,q)  ^ 

(66) 


All  the  variables  in  Eqn.  (66)  have  been  nondimensionalized  using  the  same  scales  introduced  in 
previous  traveling  wave  analysis.  In  addition  we  have  nondimensionalized  velocities  v  and  u  with 
the  constant  traveling  wave  speed,  !7,  and  the  number  density,  ti,  with  n©  =  (t>o/Vo  where  Vo  is  the 
initial  volume  of  a  bubble.  As  with  the  case  of  no  relative  motion,  there  are  three  main  parameters 
controlling  the  qualitative  structure  of  the  shock  wave.  Fig.  (5)  shows  the  solutions  to  the  ODEs 
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represented  in  Eqn.  (66)  for  a  viscous  damping  of  C  =  0.1  and  a  range  of  values  for  the  Mach 
number,  M,  and  initial  volume  fraction  of  gas,  <^o. 


M=1.01  M=1.1  M-1.5 


Figure  5: 

The  plots  in  Fig.  (5)  are  arranged  in  a  matrix  form,  with  solutions  of  the  same  Mach  number  in 
columns  and  solutions  with  the  same  initial  gas  fraction  (ahead  of  the  shock)  in  rows.  From  the  nine 
plots  shown  in  Fig.  (5)  we  can  see  three  distinct  shock  waveforms.  The  shock  in  Fig.  (5b)  begins  with 
a  sharp  rise  in  pressure  to  an  over-peaked  value  followed  by  a  oscillatory  region.  From  Fig.  (5f)  we 
again  see  a  sharp  rise  in  pressure  with  an  oscillatory  region,  but  these  oscillations  are  about  a  lower 
pressure  value  than  the  final  equilibrium  value.  The  oscillations  are  followed  by  a  relaxation  region 
which  is  due  solely  to  the  relative  motion  of  the  gas  bubbles  and  not  their  oscillations  (the  bubbles 
have  stopped  oscillating  by  this  point).  The  last  qualitative  waveform  mentioned  here  is  illustrated  in 
Fig.  (5g).  This  waveform  results  in  a  slower  monotonic  rise  in  pressure  with  no  oscillatory  behavior. 
The  suppression  of  oscillations  in  the  shock  waveform  was  a  direct  result  of  allowing  the  bubbles  to 
move  relative  to  the  liquid  in  the  model.  These  results  are  consistent  with  the  results  of  Ishii  et  al 
[9]. 
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5  Conclusion 


We  have  developed  a  non-equilibrium,  nonlinear  equation  of  state  (EOS)  that  provides  a  dynamic 
relation  between  pressure  in  the  mixture  and  the  mixture  density  and  the  number  density.  The 
EOS  contains  the  first  two  material  time  derivatives  of  both  the  mixture  density  and  the  number 
density,  allowing  for  the  possibility  for  solutions  with  oscillatory  behavior.  We  have  examined  some 
of  the  possible  traveling  wave  solutions  obtained  when  the  nonlinear  EOS  (with  and  without  relative 
motion)  is  combined  with  the  fully  nonlinear  equations  of  mass  and  momentum  conservation.  The 
structures  of  the  shocks  solutions  were  found  to  agree  qualitatively  with  the  waveforms  observed  in 
the  experiments  performed  by  Noordzij  and  van  Wijngaarden.[2] 

The  attenuation  of  pressure  waves  are  related  to  the  damping  mechanisms  that  exist  for  bubbles. 
The  three  most  important  damping  mechanisms  are  associated  with  heat  transfer  between  the  gas 
bubble  and  the  liquid,  the  drag  due  to  relative  motion  of  the  bubble  and  viscous  dissipation  at 
the  bubble  interface.  We  have  incorporated  the  latter  two  damping  mechanism  in  our  model.  In 
addition,  in  our  research  we  have  shown  how  one  can  introduce  an  effective  damping  parameter 
that  captures  the  thermal  dissipation  which  occurs  in  a  bubble  oscillating  periodically.  We  have 
shown  that  the  effective  damping  parameter  allows  a  more  accurate  description  than  that  obtained 
by  modeling  the  bubbles  as  either  purely  isothermal  or  adiabatic. 
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